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Abstract 

We design and analyze an adaptive /ip-finite element method (hp-AFEM) in dimen¬ 
sions n = 1,2. The algorithm consists of iterating two routines: hp-NEARBEST finds 
a near-best hp-approximation of the current discrete solution and data to a desired accu¬ 
racy, and REDUCE improves the discrete solution to a finer but comparable accuracy. 
The former hinges on a recent algorithm by Binev for adaptive hp-approximation, and 
acts as a coarsening step. We prove convergence and instance optimality. 


1 Introduction 

The discovery that elliptic problems with localized singularities, such as corner singu¬ 
larities, can be approximated with exponential accuracy propelled the study and use of 
hp-FEMs, starting with the seminal work of Babuska. The a priori error analysis origi¬ 
nated in the late 70’s with the earliest attempts to study the adaptive approximation of 
a univariate function, having a finite number of singularities and otherwise smooth, by 
means of piecewise polynomials of variable degree [HI Eg. These results influenced Gui 
and Babuska in their pioneering study of the convergence rate of the hp-approximation to 
a one dimensional model elliptic problem in |26j and in their subsequent work Eg, which 
proves convergence of an adaptive hp-algorithm with a predicted rate. However, due to 
the assumptions on the admissible error estimators, which appear to be overly restrictive, 
the results in [Sg cannot be considered completely satisfactory. Starting from the late 
80’s the study of a posteriori error estimators and the design of adaptive hp-algorithms 
has been the subject of an intense research. We refer to the book [3g, and the survey 
paper |15] . as well as the references therein for more details. 

However, despite the interest in hp-FEMs, the study of adaptivity is much less de¬ 
veloped than for the h-version of the FEM, for which a rather complete theory has been 
developed in the last decade Ea Eg El [mill]; we refer to the survey m- Regarding the 
hp-FEM, we mention miEiiigEi which prove convergence without rates. The purpose 
of this paper is to bridge this gap: we present a new hp-AFEM, which hinges on a recent 
algorithm by Binev for adaptive hp-approximation mg, and prove several properties 
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including instance optimality in dimensions n = 1,2. The theory is complete for n = 1 
but there are a couple of pending issues for n = 2, which we discuss below. 

The success of hp-AFEM’s hinges on having solutions and data with suitable sparsity 
structure, as well as practical algorithms that discover such a structure via computation. 
This is why existing hp-AFEM software typically probes the current discrete solution 
to learn about the local smoothness of the exact solution, but can only search around 
the current level of resolution. We refer, in particular, to the algorithms presented in 
[SI 111311 [31 El for strategies based on analyticity checks or local regularity estimation 
(see also jH (H), to the algorithms in [31 III UHl El and [3H] for strategies based on 
the use of suitable discrete reference solutions, and to the algorithm in [35] for a strategy 
based on comparing estimated and predicted errors. 

1.1 Challenges of /ip-Approximation 

To shed light on the difficulties to design hp-AFEM, we start with the much simpler 
problem of hp-approximation for n = 1. Let := (0,1) and A be a dyadic interval 
obtained from Kq = Cl. Let p be the polynomial degree associated with A at a certain 
stage of the adaptive algorithm, and denote D = {K,p). Given v G L^(Cl) and p > 0, let 

en(v) := min ||u - and Qd(v) := argmin jjv - (1.1) 

iperp(K) ifiePp(K) 

the latter function being extended with zero outside A. The following algorithm gen¬ 
erates a sequence of hp-decompositions (I>f)^o corresponding piecewise polynomial 
approximations V£ = v-Dg. With vg := Qko,o('v), for £ > 0 and any D = (K,p) G T>£, 

• compute eK,p+i(v — vi) as well as eK',p{v — vi) and ok", p(v — vi) for A' and A" 
being the two children of A; 

• if e.K,p+i(y — vi) < eK',p(v — ve) + eK",p{v — vi), then replace D hy D := {K,p -|- 1) 
in D^+i and set 

vi+i := vi + Q^{v - ve); 

• otherwise, replace D by D' := {K',p) and D” := {K”,p) in and set 

vi+i := vt + Qd'{v - vi) + Qd"{v - Vf). 

Although this algorithm is deliberately very rudimentary so as to simplify the discussion, it 
mimics existing schemes that query whether it is more advantageous to refine the element 
A or increase the polynomial degree p by a fixed amount, say 1. We wonder whether 
such an algorithm may lead to near-optimal hp-partitions. In order to elaborate on this 
question, we now present two extreme examples that illustrate the role of sparsity for the 
design of hp-AFEM. 

Example 1: Lacunary Function. For a given integer L > 0, let u be a polynomial of 
degree p := 2^ — 1, such that, on each dyadic interval A of generation 0 < i < L, v is 
L^-orthogonal to the linear polynomials with vanishing mean. Since we need to impose 2^ 
orthogonality relations for each level i, we get altogether 1 -|- 2 -|- 2^ -|- • • • -I- 2^“^ = 2-^ — 1 
constraints. We thus realize that a nontrivial polynomial of degree p does exist because 
it has 2^ parameters. We also see that the algorithm above bisects all dyadic elements 
A starting from Kg until reaching the level L, and that V£ for all 0 < £ < L is the 
piecewise constant function that takes the mean-value of v on each element in Di. Even 
if the algorithm stops refining at level L and chooses from then on to raise the polynomial 
degree by 1 in each of the p elements, then at least p new degrees of freedom have to 
be added in each interval to represent v exactly. This leads to a total of p^ degrees of 
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freedom activated to capture a polynomial of degree p, thereby proving that this process is 
non-optimal. We conclude that to be near-optimal, hp-AFEM must be able to backtrack 
and review decisions made earlier. This process, from now on called coarsening, is missing 
in most algorithms for hp-adaptivity except, for example, that of Demkowicz, Oden and 
Rachowicz [20], for which there are no optimality results. The preceding function is 
extremely sparse for hp-approximation, in fact a single polynomial, but its structure is 
hard to discover in practice because of the sparsity gap. 


Example 2: Non-degenerate Function. We now consider the canonical function 
v(x) = a;“ with a < 1 on 0 = (0,1), studied by DeVore and Scherer [22] and by Gui and 
Babuska [26], which does not exhibit a sparsity gap. In fact, the following non-degeneracy 
property is valid: there exist constants Ci, C 2 such that for all intervals K and polynomial 
degrees p 


^ / (iK,p+l(v) 

E2 S - 

(iK,p[v) 


< Cl. 


The exponential rate of convergence derived a priori in HZj, as well as the linear increase 
of polynomial degrees starting from the origin, are based on this crucial property. Similar 
results have been derived later for n = 2 by Babuska and Guo (SSJ [22] and for n = 3 
by Schotzau, Schwab and Wihler [421143| ; see |39| for a thorough discussion of the cases 
n = 1,2. It is thus conceivable, as observed in practice, that decisions made by hp- 
AFEM’s with a building block such as that above do not produce unnecessary degrees of 
freedom for problems such as Example 2. The lack of a coarsening step in most existing hp- 
software could thus be attributed to the very special geometric features of point and edge 
singularities, this being a special rather than a universal property to design an optimal 
hp-AEEM. 


1.2 Our contributions 

Since we wish to account for a large class of functions (solutions and data), perhaps ex¬ 
hibiting degeneracies such as in Example 1, our hp-AEEM includes a coarsening routine, 
which we envisage to be unavoidable for obtaining optimality. Our hp-AEEM hinges on 
two routines, hp-NEARBEST and REDUCE, and the former in turn relies on the adap¬ 
tive hp-approximation routine by Binev mis]. To describe them, let u = u{f, A) S Hq (O) 
be the solution to a second order elliptic PDE on a domain O C K",n = 1,2, with data 
(/, A), where / denotes forcing term(s) and A parameters such as coefficients. 

Given a reduction factor g G (0,1), a conforming hp-partition D, (discontinuous) 
hp-FEM approximations (/d,Ad) to (/, A) over B, the routine REDUCE produces a 
conforming hp-refinement D such that the | • -error in the (continuous) Zip-fern Galerkin 
approximation on CD to the exact solution u(/d, Ad) is less than p times the same Galerkin 
error relative to the partition B. This routine will be implemented as an AFEM routine 
that applies under a no-data-oscillation assumption. 

The routine hp-NEARBEST deals with nonconforming meshes and subordinate dis¬ 
continuous functions. Given a tolerance e > 0, a generic function v G and data 

(/, A), hp-NEARBEST produces a nonconforming hp-partition B and suitable projec¬ 
tions (/d, Ad) of the data onto discontinuous hp-FEM spaces over B. The output is such 
that the square root of a specific error functional is less than s. This error functional 
is defined as the sum of the squared broken | • |//i-error in the best (discontinuous) Zip- 
approximation over B to u and 5“^ times the squared Zip-data oscillation osc|,(/, A) over 
B, for a sufficiently small penalty parameter i5 > 0. In turn, osc|,(/. A) measures the er¬ 
rors / — /d and A — Ad on the partition B in such squared local norms, that the following 
bound, expressing the continuous dependence on data of the underlying linear problem, 
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holds: 


\u — m(/d, AD)|//i(n) < oscd{/, A). (1.2) 

The procedure hp-NEARBEST is based on Binev’s algorithm and is instance optimal 
for this functional. 

Our algorithm hp-AFEM consists of a repetition of calls of hp-NEARBEST and 
REDUCE with decreasing error tolerances. The calls of hp-NEARBEST, with V being 
the current approximation to the solution u, are made to guarantee instance optimality 
of the coarsened approximations. Coarsening, however, increases the error by a constant 
factor. This must be compensated by a judicious choice of the reduction factor g of 
REDUCE so that the concatenation of the two routines produces a converging sequence. 
To realize this idea we must account for the following additional issues. 

Making meshes /ip-conforming: After a call of hp-NEARBEST, the generally non- 
conforming hp-partition D has to be refined to a conforming one 6(1)) so that it can serve 
as input for REDUCE. This is obviously an issue for dimension n = 2 but not for n = 1, 
in which case one can take C(D) = D. One may wonder whether the cardinality of C(D) 
can be bounded uniformly by that of 2) for n = 2. To see that the answer is negative 
in general consider the following pathological situation: a large triangle of D with high 
polynomial degree is surrounded by small triangles with polynomial degree 1. This is the 
reason why, without further assumptions on the structure of the solution it, we cannot 
guarantee for n = 2 an optimal balance between the accuracy of the /ip-approximations 
and the cardinality of the /ip-partitions at stages intermediate to consecutive calls of 
hp-NEARBEST. Resorting to a discontinuous hp-AFEM would cure this gap at the 
expense of creating other difficulties. 

Making functions continuous: In order to quantify the reduction factor g of RE¬ 
DUCE we must be able to compare the (broken) i?^(D)-errors of the best continuous and 
discontinuous /ip-FEM approximations over 6(2)). We show that the former is bounded 
by the latter with a multiplicative constant which depends logarithmically on the maximal 
polynomial degree for n = 2. This extends upon a recent result of Veeser for the /i-version 
of the FEM Such constant does not depend on the polynomial degree for n = 1. This 
construction is needed for the analysis of hp-AFEM only but not its implementation. 


Dealing with a perturbed problem: When, preceding to a call of hp-NEARBEST, 
the current (continuous) hp-approximation to u has a tolerance e, hp-NEARBEST will 
be called with a tolerance ~ e in order to guarantee optimality of the coarsened discon¬ 
tinuous /ip-approximation. In addition, hp-NEARBEST produces new approximations 
(/d,Ad) to the data to be used in the subsequent call of REDUCE. The prescribed 
tolerance ensures, in view o f the definition of the error functional, that oscd(/, A) ^ y/S s. 
Hence, concatenating with (1.2), we are guaranteed that |u — ii(/d, AD)|j:/i(f2) < \/ 6 e . The 
routine REDUCE approximates the solution ii(/d,Au), and so cannot be expected to 
produce an approximation to u that is more accurate than ii(/d. Ad). Therefore, in order 
to obtain convergence of the overall iteration, the condition |it — ii(/d, AD)|iri(n) < is 
needed for some parameter ^ € [0,1), which is achieved by selecting the penalty parameter 
(5 to be sufficiently small. 


The routine REDUCE will be implemented as an AFEM consisting of the usual 
loop over SOLVE, ESTIMATE, MARK, and REFINE. For n = 1, we construct 
an estimator that is reliable and discretely efficient, uniformly in p. Consequently, the 
number of iterations to achieve some fixed error reduction g is independent on the maximal 
polynomial degree. 

For n = 2, we employ the residual-based a posteriori error estimator analyzed by Me- 
lenk and Wohlmuth [3^, which turns out to be p-sensitive. We show that in order to 
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achieve a fixed error reduction, it suffices to grow the number of iterations more than 
quadratically with respect to the maximal polynomial degree. This sub-optimal result is 
yet another reason for optimality degradation at stages intermediate between two consec¬ 
utive calls of hp-NEARBEST. Nevertheless, our result improves upon a recent one by 
Bank, Parsania and Sauter [3] , which requires the number of iterations to be proportional 
to the fifth power of the maximal polynomial degree. 

Throughout this work, we assume that the arising linear systems are solved exactly. 
To control the computational cost, optimal iterative solvers, uniformly in the polynomial 
degree would be required. We refer to [S] for an example. 

This work is organized as follows. We present hp-AFEM within an abstract setting in 
Sect, [^and prove that it converges, and that the sequence of outputs of hp-NEARBEST 
is instance optimal. We give a brief description of Binev’s algorithm in Sect. In Sect. 
1^ we apply the abstract setting to the general 1-dimensional elliptic problem. Finally, in 
Sect, [^we apply the abstract theory to the Poisson equation in two dimensions. 

The following notation will be used thoughout the paper. By 7 ^ <5 we will mean that 
7 can be bounded by a multiple of <5, independently of parameters which 7 and 6 may 
depend on. Likewise, 7 ^ ^ is defined as <5 < 7 , and 7 ~ ^ as 7 < 5 and 7 > i5. 


2 An abstract framework 

We now present the hp-AFEM in two steps. We first deal with an ideal algorithm 
and later introduce a practical scheme including REDUCE. We also discuss a possible 
realization of REDUCE. 

2.1 Definitions and assumptions 

On a domain O C ffi", we consider a, possibly, parametric PDE 

Axu = f. (2.1) 

Here the forcing / and the parameter A (representing, e.g., the coefficients of the operator) 
are taken from some spaces F and A of functions on O, such that there exists a unique 
solution u = u{f, A) living in a space V of functions on O. We assume, for simplicity, that 
V and F are Hilbert spaces over M. 

We assume that we are given an essentially disjoint initial partition of H into finitely 
many (closed) subdomains (the ‘element domains’). We assume that for each element 
domain K that we encounter, there exists a unique way in which it can be split into element 
domains K' and K", the ‘children’ of K, such that K = K' \J K" and \K' (7 K"\ = 0. The 
set .S of all these element domains is therefore naturally organized as an infinite binary 
‘master tree’, having as its roots the element domains of the initial partition of H. A finite 
subset of ^ is called a subtree of the master tree when it contains all roots and for each 
element domain in the subset both its parent and its sibling are in the subset. The leaves 
of a subtree form an essentially disjoint partition of H. The set of all such ‘/i-partitioi^’ 
will be denoted as K. For IK, 3C € K, we call % a refinement of X, and denoted as 3C < IK, 
when any AT S IK is either in X or has an ancestor in X. 

Our aim is to compute ‘/ip-finite element’ approximations to u, i.e., piecewise polyno¬ 
mial approximations, with variable degrees, w.r.t. partitions from K. In order to do so, it 
will be needed first to replace the data (/, A) by approximations from finite dimensional 
spaces. For that goal as well, we will employ spaces of piecewise polynomials, with variable 
degrees, w.r.t. partitions from K, as will be described next. 
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For all K G let Vk, Fx, be (infiiiite dimensional) spaces of functions on K, such 
that for any DC G K, it holds that, possibly up to isomorphisms, 

VC YIvk, f= H Fk, a C J] Ak C a. 

ifeJC ifeac kgx 

Here A is a subset of A, which contains all the parameters that will be allowed in our 
adaptive algorithm hp-AFEM, and, for simplicity, has a Hilbert topology. For all {K, d) G 
X N (hereafter N stands for the set of all positive natural numbers) and Z G {V, F, A}, we 
assume finite dimensional spaces Zx,d C Zx of functions on K such that Zx,d V Zx^d+l^ 
Zx,d C ZxGd ^ Zx"^d-/ and Z n ^K^d Is dense in Z . 

In applications, Vx,d will be a space of polynomials of dimension ~ d. For instance, 
when K is an n-simplex, Vx,d may be chosen as Vp{K), where the associated polyno¬ 
mial degree p = p{d) can be defined as the largest value in N such that dimPp_i(iC) = 
F: d. This definition normalizes the starting value p(l) = 1 for all n G N. Only 
for n = 1, it holds that p{d) = d for all n G N. 

Analogously, the spaces Fx,d and Ax,d will be selected as (Cartesian products of) 
spaces of polynomials as well, of degrees equal to p plus some constant in Z. 

In the following, D G ^xN will denote an hp-element: it is a pair {Kx, do) consisting 
of an element domain Kx G A, and an integer G N. We will write Zx = Zxo,dD- 

For all Z? G .ft X N, we assume a projector Qx : V x F x A ^ Vx x Fx x Ax, 
and a local error functional ex = ex{v, /, A) > 0, that, for {v,f,X)GV x F x A, gives a 
measure for the (squared) distance between {v\xo: f\Kn^ Mkh) local approximation 

(vx, Id, Ad) := Qx{v, f, A). We assume that this error functional is non-increasing under 
both ‘h-refinements’ and ‘p-enrichments’, in the sense that 

ex' + ex" < ex when Kx', Kx" are the children oi Kx, and dx' = dx" = dx', (^2 2) 
cd' < cd when Kx' = Kx and dx' > dx- 

A collection CD = {D = {Kx,dx)} of /ip-elements is called an hp-partition provided 
DC(®) := {Kx '■ D G ®} G K. The collection of all /ip-partitions is denoted as D. For 
2) G D, we set the /ip-approximation spaces 

Zx-.= Y[Zx, {Zg{V,F,A}), 

DeD 

and define 

# 2 ) := ^ dx- 
DeD 

In our applications, the quantity ^2) is proportional to the dimension of Zx, and 
ex(v, fx,Xx) is the sum of the squared best approximation error of v\ko from Vx in 
I • \m[KD) 5“^ times the square of the local data oscillation. 

For 2) G D, we set the global error functional 

Vxiv,f,X) := ^ ex{v,f,X), 

DeD 

which is a measure for the (squared) distance between {v, f. A) and its projection 

( n n fo, Ad^ G Vx X Fx X Ax- (2.3) 

Dec DeD r>eD 

For 2), CD G D, we call D a refinement of 2), and write D < D, when both 2<C(2)) < CK(CD), 
and djj > dx, for any D G D, D G D with Kx being either equal to Kp^ or an ancestor 
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of Kjj. With this notation, (2.2) is equivalent to 


E5(n,/,A) <Ea5(i;,/,A) VD > 2). (2.4) 

We will apply a finite element solver that generally operates on a subset D'^ of the set of 
hp-partitions D, typically involving a restriction to those CD G D for which the ‘h-partition’ 
3C(CD) is ‘conforming’. We assume that there exists a mapping C : D —>■ O'” such that 

6(2)) >2) V23 G D. (2.5) 


We emphasize that even for CD G D'’, generally the space Ed is not a subspace of V. 
Conforming subspaces, used e.g. in Galerkin approximations, are defined as 


VCg := Ed n E 


( 2 . 6 ) 


With regard to (2.3), we introduce the notation 


Id '■= fo, Ad := Ad, 

DeD dgd 

but reserve the symbol ud to denote later a suitable near-best approximation to u G E 
from Eg. 


2.2 A basic /ip-adaptive finite element method 

Our aim is for given (/, A) G F x A and £ > 0, to find D with an essentially minimal 
such that Ed(m(/, A),/, A) < e. We will achieve this by alternately improving either the 
efficiency or the accuracy of the approximation. To that end, we begin by considering a 
basic algorithm, which highlights the essential ingredients of a /ip-adaptive procedure. We 
make use of the two routines described below. The first routine is available and will be 
discussed in Sect. Since we are not concerned with complexity now, existence of the 
second routine is a simple consequence of the density of the union of the /ip-approximation 
spaces in E. 

• [®,/d,Ad] := hp-NEARBEST(£,u,/,A) 

The routine hp-NEARBEST takes as input £ > 0, and (u, /, A) G E x F x A, 
and outputs 2) G D as well as (/d,Ad) such that Ed(u, /, A )2 < e and, for some 
constants 0 < & < 1 < B, #2) < i?#2) for any 2) G D with Eg(u, /, A)5 < be. 

• [2i>, u] := PDE(£, 2), /d. Ad) 

The routine PDE takes as input £ > 0, 2) G D'’, and data (/d, Ad) G Fd x Ad. It 
outputs 2) G D'^ with 2) < 2) and u G Eg such that ||u(/d, Ad) — u\\v < £• 

The input argument v of hp-NEARBEST will be the current approximation to u(/. A). 
In an ‘/i-adaptive’ setting, usually the application of such a routine is referred to as 
‘coarsening’. Since the data (/d,Ad) G Fd x Ad of PDE is discrete, it will be said to 
satisfy a no-data-oscillation assumption w.r.t. 2 ). 

We make the following abstract assumptions concerning the relation between the error 
functional, the norm on E, the mapping (/, A) i—>■ u(/. A), and the constant b of hp- 
NEARBEST. We assume the existence of constants Ci,C 2 > 0 with 

C 1 C 2 < b, (2.7) 
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such that 


l|w(/,A)-u(/i,,AD)||y <Ci inf Ed(u;,/,A)5 VD G D, V(/, A) G E x A, (2.8) 

wGV 

sup f, X)^—E,-x){v, f, X)^\ < C 2 \\w — v\\v VD G D, Vu, ic G F. (2.9) 

if,\)eFxA 


The condition (2.91 means that Ed(w,/, A)! is Lipschitz w.r.t. its first argument. In 
our applications, we will verify this condition with (72 = 1. The condition (2.8) will be a 
consequence of the continuous dependence <111 of the solution on the data, and the fact 
that the error functional will contain the square of a data oscillation oscd(/. A). Since 
this term is penalized by a factor we will be able to ensure (2.8) with Ci ~ which 
yields (2.7) by taking S sufficiently small. 

Our basic hp-adaptive finite element routine reads as follows. 

hp-AFEM(uo, /, A, eo) 

% Input: (uo, /, A) G E X F X A, eo > 0 with |iM(/, A) — Mo||y < eo- 
% Parameters: > 0 such that C 1 C 2 < b{l — p), w G (^, M G (0,1), 

for i = 1,2,... do 

[Di, fDi,X'Di] :=hp-NEARBEST(a;ei_i, Mi_i, /, A) 

[Vi, Ui] := PDE(/xei_i, C{Vi), /d;, AdJ 

•= (M + CiUj)Si-i 

end do 

Note that buj — C 2 > 0, and that e^ = (/r + (7ia;)®eo, where p + CiUJ < 1. 


Theorem 2.1. Assuming (2.7)-(2.9), for the sequences {ui), (Vi) produced mhp-AFEM, 
writing u = u(f. A), it holds that 


and 


\u-Ui\\v<Si Vi > 0, Ei,,(u,/,A )2 < (w + ( 72 )ej_i Vi > I, (2-10) 


ffVi < BffV for any V G J} with Fjd(u, f, X)^ < (bu! — C 2 )£i-i- (2-11) 


Proof. The bound ||u — uo||y < eo is valid by assumption. For i > 1, the tolerances used 
for hp-NEARBEST and PDE, together with (2.8), show that 


( 2 . 12 ) 


IIU - Willy < ||w(/i,-, AdJ - Wj||y + ||m - u(f'D.,Xq:).)\\v 

^ h-Si-i + Cl E®. (ui-i, /, A) 2 < (p + Ciuj)ei-i = et. 

The first statement follows for all i > 0. Using this and ( |2.9[ ) implies the second assertion 
Ed,(u,/,A) 5 < Ei,,(itj_i,/, A)5 + C' 2 ||w - Ui-i||y < (uj + C 2 )ei-i Vi > 1. 


Let 2) G D with Ed (it,/, A) 2 < (ftw — ( 72 )ei-i- Then, again by (2.9), ED(wi-i,/, A)2 < 
buJEi-i and so ffVi < BffV because of the optimality property of hp-NEARBEST. □ 


The main result of Theorem |2.1| can be summarized by saying that hp-AFEM is 
instance optimal for reducing ED(it(/, A),/, A) over 2) G D. Recall that in our applica¬ 
tions, Ed(u(/, A), /, A) will be the sum of the squared best approximation error in u from 
the nonconforming space Ed = YIdeD in the broken TJ^-norm and a squared data 
oscillation term penalized with a factor <5“^. 

Additionally, Theorem 2.1 shows linear convergence to u of the sequence (wj) of 
conforming approximations, in particular Ui G Eg. where V^ B Vi > C(Vi). Since 












Ei = [^ + C'iw)*eo, the infinite loop in hp-AFEM can be stopped to meet any desired 
tolerance. 

The preceding algorithm hp-AFEM has the minimal structure for convergence and 
optimality. Since the routine PDE neither exploits the current iterate nor work already 
done, we present a practical hp-AFEM in Sect. 2.3 which replaces PDE by REDUCE. 


Finally in this subsection, we comment on the implications of the instance optimality 
result concerning class optimality. For A € N, let := argmin{Eii(u,/, A)^ : D € 
D, #2) < A} and let the best approximation error be 


(Tat := (u,/, A) =. 


Remark 2.1 (algebraic decay). If crjv decays algebraically with A, namely sup^r N^aN < 
oo, then for the sequence (Vi) produced in hp-AFEM, one infers that Ed.(u,/, A)5 
decays algebraically with #1)^ with the optimal rate: supj(#Di)® Ed-(m,/, A) s < oo. In 
other words, instance optimality implies algebraic class optimality. 

Remark 2.2 (exponential decay). For hp-approximation, it is more relevant to consider 
an exponential decay of ctat, i.e., sup at < oo for some ri,T > 0. This is precisely 

the situation considered in [niiiaiia] for adaptive Fourier or Legendre methods. 

Let us asssume, for convenience, that = C^e~'^^ for some constant and 
ignore in subsequent calculations that A has to be an integer. In view of Theorem |2.1[ 
let A and e^-i be so related that = {buj — C 2 )£i-i Since apparently #1)^ < BN and 
Ed-(u,/, A)5 < (w + C' 2 )ei-i, we deduce 

8 up(e’i»^-)'E„,(„,/,A)^) 


with fi := B~'^ri. 

On the other hand, we will see in Corollary |3.1| that the routine hp-NEARBEST 
satisfies its optimality conditions for any R > 1, at the expense of 6 = b{B) 0 when 
R 1. Moreover, as we have seen, in our applications we will be able to satisfy (2.7)- 
(2.9) for any 6 > 0 by taking the penalization parameter 5 small enough. Therefore, we 
conclude that if ctat decays exponentially, characterised by parameters {r],T), then so do 
the errors produced by hp-AFEM for parameters where fj = B ~'^ r ] can be chosen 

arbitrarily close to ry (at the expense of increasing the supremum value). This situation is 
much better than that encountered in 


2.3 The practical /ip-adaptive finite element method 

To render hp-AFEM more practical we replace the routine PDE by REDUCE, which 
exploits the work already carried out within hp-AFEM and reads 

• [D,u] :=REDUCE(e,'D,/D,AD) 

The routine REDUCE takes as input a partition T) G D'^, data (/d, Ad) S Rd x Ad, 
and a desired error reduction factor g G (0,1], and produces a conforming partition 
2) = 2)(23, g) G D'^ with 2) > 2) and a function u G Ug such that 

||■w(/D, Ad) - hlly < £» inf ||u(/d, Ad) --cHv- (2.13) 

vev^ 

Inside the practical hp-AFEM, the routine REDUCE will be called with as input 
partition the result of mapping C : D —>■ applied to the output partition of the preceding 
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call of hp-NEARBEST. In order to bound the right-hand side of (2.13), we make the 
following assumption: 


inf Wv — wWv^CzT) inf V® G D, Vi; S E. (2-14) 

’ (/,A)eFxA 

In our applications, the infimum on the right-hand side reads as the squared error in the 
broken R^-norm of the best approximation to v from = Wd^d^d- The left-hand 
side reads as the squared error in of the best approximation to v from = 

Rg(r2)nn£)gD ^D- The constant should ideally be independent of ®. We will see in 
Sect. 1^ that this is the case for our application in dimension n = 1. However, for n = 2 we 
will show in Sect. that depends logarithmically on the largest polynomial degree; 
this extends a result by A. Veeser [15] . 

Our practical hp-adaptive finite element routine reads as follows: 

hp-AFEM(Mo, /, A, eo) 

% Input: (uo, /, A) G E X F X A, Eg > 0 with ||m(/. A) — Molly < eo- 
% Parameters: p, w > 0 such that C 1 C 2 < b{l — p), w G (^, p G (0,1). 
for i = 1, 2,... do 

[Vi, /cDi, AdJ :=hp-NEARBEST(a;ei_i, Mi_i, /, A) 

[®„ M,] := REDUCE(^^(p^J^j^^, e(®.), /d., AdJ 

Si := (p -|- Ciuj)ei-i 

end do 


Corollary 2.1 (convergence and instance optimality). Assuming (2.7)-(2.9) and (2.14), 
the sequences (ui), (®i) produced in the practical hp-AFEM above satisfy properties 


(2.10) and (2.11) in Theorem 2.1 


Proof. In view of the second part of the proof of Theorem 2.1 it is sufficient to prove that 
||u ~ fti||y < e^. We argue by induction. If ||m — Mi||v < e^-i, which is valid for i = 1, 
then, after the ith call of hp-NEARBEST, (2.14) and (2.8) imply that 


inf ||m(/d,,Ai5 ,) -M||y 

< ||u - Mi_i||y -I- inf ||Mi_i -M||y-I-|| u-m(/d,,AdJ||v 

< Si -1 + ^DiiUi-l, f, A) -I- Cl Ex>,(ui- 1 , /, A) 

— (1 + (C'a.Di + Ci)uj)ei-i. 


(2.15) 


Consequently, after the subsequent call of REDUCE, it holds that ||m(/i 3 , , AdJ — Mi||y < 
HEi-i according to (2.13). This result combined with (2.12) shows that ||m —Mi||v < Si. □ 


Remark 2.3 (complexity of hp-AFEM). Let us consider the case that the constants 
Ca.D, defined in (2.14), are insensitive to T), namely, 


C 3 := sup < 00 . 
dgd 


(2.16) 


This entails that the reduction factor Qi = ^ of REDUCE satisfies inf^ Qi > 

0. Additionally, suppose that, given a fixed g G (0,1], REDUCE realizes (2.13) with 


sup 


#T>(^,g) 

#® 


< 00 . 


(2.17) 
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If, furthermore, 


(2.18) 


r, #e(®) / 

^4 := sup < oo, 

dgd 

then the sequences (CDi)^ and {T>i)i produced in hp-AFEM are so that In 


view of the optimal control over #I)i, given by Theorem 2.1 and Corollary |2.1[ we would 
have optimal control over the dimension of any /ip-finite element space created within 
hp-AFEM. This ideal situation only happens in the one-dimensional case. 


2.4 A possible realization of REDUCE 

Let A\ G £(V,V') for all A S A and define the associated continuous bilinear form 
a\{v,w) := {A\v,w) for any v,w G V, where (•,•) denotes the duality pairing between 
V and V. We assume that A\ is symmetric, which is equivalent to the symmetry of 
the form a\. We furtherly assume that each a\ is continuous and coercive on V, with 
continuity and coercivity constants a* > a* > 0 independent of A G A. It is convenient 
to introduce in V the energy norm |||u|||a = a/ ax{v, v) associated with the form ax, which 
satisfies ycU||?^||y < |||r^|||A < ■\/A*||?^||y for all v G V. Let F c V. 

Given D G D and data (/, X) G F x A, the (Galerkin) solution md(/. A) G Idg of 

axMf,X),v) = {f,v) Vug Kg (2.19) 

is the best approximation to u{f,X) from Kg in ||| • |||a. In view of a posteriori error 
estimation, we will consider Galerkin solutions from Kg only for data in Eg x Ad, i.e., 
for data without data oscillation w.r.t. D. 

For D G D°, D G D, let us introduce local a posteriori error indicators 


■ Kg X Ed X Ad —>■ [0,oo), 


which give rise to the global estimator 


£-d{v, /d, Ad) 



( 2 . 20 ) 


Given data (/d)Ad) without data oscillation w.r.t. D, £d(' 1 ' 5 /dj Ad) will be used with 
V = ud(/d!Ad) as an estimator for the squared error in this Galerkin approximation 
to u(/d,Ad)- It should not be confused with Ed(u, /, A), the latter being the sum of 
local error functionals eoiv, f, X), that estimates the squared error in a projection on 
Kd X Ed X Ad of (u, /, A) G K x E x A. 

Given any M C 2), it will be useful to associate the estimator restricted to M 

£d(M,u,/d, Ad) := f X/ 

Vdgm / 


We assume that £d satisfies the following assumptions: 

• Reliability: For D G D'^, and {fs, Ad) e Ed X Ad, there holds 


||'a(/D, Ad) - ud(/d, AD)||y < £d(ud(/d, Ad), /d. Ad)- ( 2 - 21 ) 


• Discrete efficiency: For D G D'^, (/d. Ad) G Ed x Ad, and for any M C T>, there 
exists a 2)(M) G with !f)(M) > 2) and < #23, such that 

II^d(m)(/d,-^b) - wb(/b, AD)||y > £©(^1 ,ud(/d. Ad),/d. Ad)- (2.22) 
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Then a valid procedure REDUCE is defined as follows. 


[i),U£] = REDUCE(p,D,/i,,Ai,) 

% Input: g G (0,1], T) G D°, (/d, Ad) G x A®. 

% Output: T) gW with D > D, and the Galerkin solution Uf , = md(/d, Ad). 
% Parameters: 9 G (0,1] fixed. 

Compute M := M{g) gN sufficiently large, cf. Proposition 
Do := D; SOLVE: compute udo(/d,Ad) 
for i= 1 to M do 

ESTIMATE: compute (ud,.! (/d, Ad),/d, Ad) : Dg Dj_i} 

MARK: select C Di-i with 


2.1 


Ad),/d, Ad) > ^£Di_i(^Di-i(/23, Ad),/d, Ad) 

REFINE: D, := ©(M^-i) 

SOLVE: compute MDi(/D,AD) 
end 

£> := Dm; Uf, = udm(/d, Ad) 


Proposition 2.1. Assuming ( |2.21 1 and (2.22), the number M = M(g) of iterations that 
are required so that [D, ud(/d, Ad)] = REDUCE(p, D,/d, Ad) satisfies 

I|m(/d, Ad) - wd(/i>> AD)||y < e inf ||u(/d, Ad) -'cjly 

vev% 

is at most proportional to logp“^, and ffil) < #D, both independent of T> G D'^, and 
(/d,Ad) G Fd X Ad- So both (2.13) and (2.17) are realized. 

Proof. Since /d and Ad are fixed, for simplicity we drop them from our notations. Ap¬ 
plying (2.22) with D = Di_i and D^ = D(Mi_i), the definition of MARK, and (2.21) 
we get 

IjuDi — Ul)i_i\\v ^ £Di_i(^i-l,^I>i-i>/l>> Ad) 

> /d, Ad) 

>6»||u-UD,_J|y. 

This and the uniform equivalence of || ■ ||y and ||| • |||aj, =: ||| ■ ||| give the saturation property 


IjuDi - UD,_1 III > - MDi_i III 

for some positive constant C*. Then, using Pythagoras’ identity 

|||m - WDi IIP = |||m - UDi_i IIP - llluDi - MDi-i ||P, 

we obtain the contraction property 

|||m — UDi III < «:|||u — MDi_i III 

for K = >/l — (7*0 < 1. We conclude that 


(2.23) 

(2.24) 

(2.25) 


\u - un^Wv < 


1 


1 


|||m- mdmIII < 


1 


:K^|||U - UdII 


inf |||u —u|||<i/—inf \\u — v 


o’" 




vev- 


a* 


«GV.S 


V- 
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Enforcing < Q yields M = O(loggi ^). In addition, since #Di < #®i-i 

1 <i< M according to (2.22), the proof is complete. 


for 

□ 


Remark 2.4. The partition ©(M) can be built by an ‘h-refinement’ or a ‘p-enrichment’, 
or both, of the elements D G M, if necessary followed by a ‘completion step’ by an 
application of the mapping C in order to land in O'”. The estimate #I)(M) < shows no 
benefit in taking 0 < 1, i.e., in taking a local, ‘adaptive’ refinement. In our algorithm hp- 
AFEM, the adaptive selection of suitable hp partitions takes place in hp-NEARBEST. 
Nevertheless, in a quantitative sense it can be beneficial to incorporate adaptivity in 
REDUCE as well, by selecting, for a 0 < 1, a (near) minimal set M C Di-i in MARK. 


Remark 2.5. The discrete efficiency of the estimator implies its “continuous” efficiency. 
Indeed, taking M = D in (2.22) and denoting (D = (D(D), and temporarily dropping /d 
and Ad from our notations, we have 


£-d{ud)^ < a*||uD - udWy ^ lll^h “ '“dIIIa = lll« “ '“dIIIa “ 11^^ “ '“hlllA < 11^^ “ '“dIIIa 


= inf |||u - u||U < a* mf ||m - v\\y. 


Consequently, recalling (2.21), a stopping criterium for REDUCE could be defined as 
follows 

£iii(MDi(/D, Ad),/d, Ad) < Cp£d(ud(/d, Ad),/d, Ad), 


where C is a constant in terms of the “hidden constants” in ( 2 . 21 ) and ( 2 . 22 ), and a* and 
a*. 


Assumptions (2.21)-(2.22) about reliability and discrete efficiency can be substituted 
by the following three assumptions concerning the estimator. This will be used for our 
application in two dimensions in Sect. 

• Reliability and efficiency: For D G ID)“, there exists R-rijr-s > 0, such that for 
(/d,Ad) G Ed X Ad, and ||| • Iai, =: ||| • ||| one has 


I’D £-%{uTi{fT), /d. Ad) < |||u(/d. Ad) — wd(/d, Ad)|||^ 


Stability: For D G 


< i?D £d(i^I'(/d, Ad), /d, Ad); 
and all (/d. Ad) G Ed x Ad, v,w G Vf, one has 


a/Ed £d(v, /d. Ad) — /d. Ad) 


< u — w 


(2.26) 


(2.27) 


• Estimator reduction upon refinement: There exists a constant 7 < 1, such 
that for any M C fD G D'”, there exists a T)(M) G D'^ with !f)(M) > D, #'D(M) < 
such that with § := {D G ©(M) : BD G M with Kp, C Ad}, 


Ad),/d, Ad) < 7£d(^i^d),/d, 

\ wd(/d. Ad), /d. Ad) < £d(B \ M, ud(/d. Ad), /d, Ad), 

for any (/d,Ad) G Ed x Ad. 


(2.28) 


With 0 from REDUCE and 7 from (2.28), we set 7 := (1 — 0) + Oj. For B < D G D'”, 
and (/d, Ad) G Ed x Ad, we define the squared total error to be 


■= £|i(u5(/D, Ad), /d. Ad). 


D '^D'' D' 
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Proposition 2.2. Assume (2.26), ( |2.27| ), and (2.28), and, msirfe REDUCE, takeT){3vV) 
as defined in ( 2.28| ). Let D € D'^, and (/d,Ai)) £ E® x Then consecutive iterands 
produced in REDUCE(gi, 2 ),/d, Ad) satisfy 


£D,(MDi(/D, Ad),/d, Ad) < 


1 - 


{l-Vir rv. 


2 RDi-i- 

Furthermore, for 2) G and (/d, Ad) G Ed x Ad, 


'SDi_i(MDi_i(/D, Ad), /d, Ad). 


|||u(/d, Ad) — ud(/d, Ad)||P < 2;d(^I'(/i>, Ad),/d, Ad) < 2|||m(/d, Ad) — '«d(/d, Ad)|P- 


IS 


Therefore, if sup Ed < oo and inf td > 0, then the statement of Proposition \2.1 
DgD^ DgD” 

again valid. 

Proof. Since both /d and Ad are fixed, we again drop them from our notations. Applying 

£D,(«D.-i)<7£D._,(«a3.-J- (2-29) 


MARK and (2.28) yields 


By virtue of (2.27), Young’s inequality, and (2.29), we have that for any C > 0, 


^ (1 + C) £Di('*^23i-i) + (1 + C ^)''’Dlil'*^2)i ~ WDi.iIlP 
< (1 + C)7 £Di_i('“I>i_i) + (1 + c ^)’"Di III^Di ~ MDi_i||P- 

By multiplying this inequality by , substituting ^ = 7“5 — 1 , and adding to 

Pythagoras’ identity (12.241), we obtain 


||M-MDj||^ + (l-V7)?'Di £Di(®i,MDj < + y/fi{l-y/j)rD, 


We resort to (2.26) to bound the right-hand side as follows in terms of an arbitrary 

fi G [0, 1] 

/3|||u- ■UD,_i||P + ((1 - fi) + V^)(l - £D,_i(®*-l,'«Di_i)- 

We now observe that the following function of p attains its minimum at 


max 


{^,((1-/3)- 


R-Oi-i 


+ 


(1 - V7)rD 

The proof of the first statement follows from 


Vi) } > /3* : 


= 1 - 


1-^7 


1 + 






(1-v^)r25 


> ■ The second 

-1 — 


statement is a direct consequence of (2.26), and the final statement follows directly from 
the first two. 

□ 


3 The module hp-NEARBEST 

In this section we describe briefly the algorithm and theory recently developed by P. Binev 
for hp-adaptive tree approximation [^, which constitutes the building block behind the 
module hp-NEARBEST. 
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3.1 /i-Adaptive Tree Approximation 


We first review the algorithm designed and studied by Binev and DeVore [7] for /i-adaptive 
tree approximation. Since, in this subsection, the local approximation spaces do not 
depend on d, temporarily we identify an element D with the element domain K]j, and D 
with the h-partition 3C(CD), the latter being an element of K. 

Recall that for any IK G K, the set of all RT G IK together with their ancestors form a 
tree T, being a subtree of the master tree Conversely, given such a subtree T, the set 
£(T) of its leaves is a partition in K. 

For the moment, we will assume that the master tree ^ has only one root. In the next 
subsection, in Remark |3.11 we will deal with the case that it has possibly multiple roots. 

For any K G A, let Ck > 0 be some local h-error functional. That means that it 


satisfies the key property (2.2), that in this /i-element setting reduces to subadditivity: 


ex' + gk" < dK 


where K' and K" denote the children of K. The corresponding global h-error functional 
reads 

Ex = ^ V3C G K. 

K^OC 

The notion of a best h-partition w.r.t. this error functional is now apparent: for iV G N, 
let 

aN := inf Ex. 

#X<Af 

This quantity gives the smallest error achievable with /i-partitions 3C with cardinality 
< N. In spite of the inf being a min, because the minimization is over a finite set, 
computing a tree that realizes the min has exponential complexity. 

A fundamental, but rather surprising, result of Binev and DeVore shows that a near¬ 
best h-adaptive tree is computable with linear complexity. A key ingredient is a modified 
local h-error functional ck defined as follows for all K € 


• ex '■= dK if K is the root; 

• J- := — H —— where K* is the parent of K and ex 0; otherwise ex = 0. 

This harmonic mean has the following essential properties: if ex dx*, then ex ~ ex, 
whereas if ex ~ dx*, then ex ~ This means that ex penalizes the lack of success 
in reducing the error from K* to AT up to a factor ^, provided ex = ex *, and always 
< 1 . 

The practical method consists of applying a greedy algorithm based on {exjxeoc- given 
an Ii-partition 3Cx, with ffOCx = N, construct %x+i by bisecting an element domain 
K G % with largest ex- It is worth stressing that if lack of error reduction persists, then 
the modified error functional ex diminishes exponentially and forces the greedy algorithm 
to start refining somewhere else. 

For ex being the squared L^-error in the best polynomial approximation on K of 
a function v, this may happen when v has local but strong singularity. The simple, 
but astute idea to operate on the modified error functionals is responsible alone for the 
following key result. 

Theorem 3.1 (instance optimality of /i-trees (3). Let the master tree A have one sin¬ 
gle root. The sequence of h-partitions {Xx)Neti given by the greedy algorithm based on 
(eif)_R'GX provides near-best h-adaptive tree approximations in the sense that 


Ex„ < 


N 

-(T, 

N-n + l 


Vn < N. 


The complexity for obtaining Xx is 0{N). 
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We can interpret Theorem 
whence N — n + 1 > N/ 2 and 


3.1 


as follows: given N let n 
Eacjv < 2cr|-N^. 


1"^] be the ceiling of N/2, 
(3.1) 


3.2 /tp-Adaptive Tree Approximation 


In this subsection, we return to /ip-approximations. An element is a pair {K, d) = 
{Ku, do), with K being the element domain, and d an integer. The local error functional 
e_D > 0 is required to satisfy (2.2), i.e., eK',d + eK",d < eK,d when K',K'' are the children 


of K, and eK,d' < when d' > d. The corresponding global hp-eiroi functional reads 
as 

Ed = y. CD VD e D. 


For G N, we set 


DeD 


ctat := inf Ed 
#Ti<N 


where #D = <^d- 

In our applications, drj is proportional to the dimension of the polynomial approxi¬ 
mation space that is applied on Kd so that is proportional to the dimension of the 
global hp-Umte element space. More precisely, given d, we take p = p{d) as the largest 
integer for which 

dimPp_i(Ar) = ^ < rf, (3.2) 

and corresponding to D = (K,d), we choose Pp((i)(iF) as approximation space. Conse¬ 
quently, for n > I, eK,d+i = ^K,d whenever p{d -|- I) = p{d). 

We describe an algorithm, designed by Binev HIS], that finds a near-best hp-partition. 


It builds two trees: a ghost h-tree T, similar to that in Sect. 3.1 but with degree dependent 
error and modified error functionals, and a subordinate hp-tree tP. The second tree is 
obtained by trimming the first one and increasing d as described in the sequel. 

Let IK G K, and let T denote its corresponding tree. For any AT G T, we denote by 
T(Ar) the subtree of T emanating from K, and let d{K, T) be the number of leaves of 7{K), 
i.e. 

d(iF,T) = #£(T(iF)). (3.3) 


The tree-dependent local hp-error functionals eif(T) are defined recursively starting 
from the leaves and proceeding upwards as follows: 

• eKifP) := cks provided K G £(T), 

• dKifJ) ■■= min{eif/(T) -f eK"{T'),eK 4 (K, 7 )} otherwise, 

where K',K” G 7 are the children of K. This local functional carries the information 
whether it is preferable to enrich the space (increase d) or refine the element (decrease 
h) to reduce the current error in K. The subordinate hp-tvee CP is obtained from T by 
eliminating the subtree 7{K) of a node AT G T whenever 


ex(Cr) = eK,d(K,r)- 
This procedure is depicted in Figure 

The hp-tree CP gives rise to an /ip-partition CD, namely the collection of /ip-elements 
D = {K,d) with K a leaf of CP and d = d{K,7). We have that ffTi = ffX, and D 
minimizes over all CD G D with 0C(D) < OC and d^ < d{KD,7) for all D G T), whence 

#5 < # 00 . 
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Figure 1: Ghost h-tiee T (left) with 10 leaves {^L{7) = 10). The label of each node K is d{K,7). 
Subordinate hp-tiee T (right) resulting from T upon trimming 3 subtrees and raising the values of d 
of the interior nodes of T, now leaves of IP, from 1 to 2, 3, and 2 respectively. 


This describes the trimming of the h-tvee T, but not how to increase the total cardi¬ 
nality of “J. To grow T, P. Binev uses a modified local hp-errov functional and a greedy 
algorithm that selects the leaf of T that would lead to the largest reduction of the hp- 
error in IP. We refer to [5] for the construction of the full algorithm for /ip-adaptive 
approximation. 

Theorem 3.2 (instance optimality of hp-tree [4l[5]). Let the master tree ^ have one single 
root. For all N G N, the algorithm sketched above constructs an hp-tree CPjv subordinate 
to a ghost h-tree T^v such that the resulting hp-partition Djv has cardinality = N 

and global hp-error functional 




2N 


N-n + 1 




Vn < N. 


In addition, the cost of the algorithm for obtaining Dn is bounded by 0( , 

and varies from 0{N log N) for well balanced trees to 0{N'^) for highly unbalanced trees. 


Binev’s algori thm gives a routine hp-NEARBEST that satisfies the assumptions 
made in Subsect. 2.2 for any B > 1 and b = — g): 

Corollary 3.1. Let B > 1. Given e > 0, let T> be the first partition in the sequence 

^ . . 1 

produced by Binev’s algorithm for which E, 


D 


< s. Then jfT) < B m.va.{ffT): I) € D, < 


^( 1 -^) 4 - 


Proof. Let D = 2)at, i.e., D is the iVth partition in the sequence, and ffD = N. For N = 1 

the statement is true, so let A > 1. Suppose there exists a D S D with < \J\{l. — j^) s 

and N > BffT). Then, with n := #i), we have < Af-f-Tn+i ^i) ^ 

Af-r^i+i ^(1 ~ From ~ s) — being a consequence of A > Bn and 

1 

B > 1, we get a contradiction with T) being the first one with Ef, < e. □ 
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Remark 3.1. In order to deal with the case that the master tree ^ has R > 1 roots, the 
following approach can be followed. 

We unify the R roots pairwise creating new element domains, each one being the union 
of two roots. When R > 2, this process has to be repeated until only one element domain 
remains, which will the new, single root. Obviously, this applies only when R is a power 
of 2. In the other case, we have to introduce at most [log 2 R] — 1 (empty) virtual element 
domains (and, formally, infinite binary trees of virtual element domains rooted at them). 
We denote the extended, single rooted master tree by 

Next, we extend the definition of eK,d as follows. At first we give a meaning to exp for 
each element domain K € A. Typically, for d G N, ex.d has the meaning of the squared 
error in the approximation of a quantity from a space of dimension d. Then a natural 
definition of exp is that of the squared error in the zero approximation. 

Considering now the elements in A.\ A, i.e., the newly created element domains, we 
distinguish between virtual and non-virtual element domains. For each virtual element 
domain, we set ex,d '■= 0 for any d G N U {0}. Finally, for each newly created non-virtual 
element domain K, being the union of K' and K" (one of them possibly being a virtual 
element domain), for d G N U {0} recursively we define 


e-x,d '■= 


min 

{d'.d"GNu{0}: d'+d"<d} 


ex',d’ 


+ eK",d"- 


Note that in the minimum at the right hand side d' or d" can or has to be zero. In that 
case, ex'.d' + e.x".d" has the interpretation of the squared error in an approximation on 
K that is zero on K' or K". 

It is easily checked that the error functional ex,d for {K,d) G .^x N satisfies (2.2), and 
Theorem\3.^ and Corollary\3.1\ apply. 


We close the discussion of the module hp-NEARBEST with the observation that 
in dimensions n > 1, Binev’s algorithm produces dp-partitions that are generally non- 
conforming. Since conformity is required by the module REDUCE, a post-processing 
step which makes the output partition conforming is required. The implementation of 
such a procedure in dimension 2 , and the analysis of its complexity, will be discussed in 
Sect. 15.11 


4 A self-adjoint elliptic problem in ID 

In this section we apply the abstract framework introduced in Sect, [^to a one-dimensional 
self-adjoint elliptic problem. 

4.1 The continuous problem and its hp discretization 

Let n := (0,1). Given /i ,/2 G L'^{Q) and iy,a G satisfying 

0 <i>^<v<v*<oo and 0 < cr < ct* < oo 

for some constants , v* and cr*, we consider the following model elliptic problem 

— {vu')' + cru = fi f 2 in , 
u(0) = m(1) = 0 , 

which can be written as in ( |2.1[ ) setting A = (z/, a), f = fi + f 2 & and 

Axu := —{izu'Y au G £(Rg (fl), 


(4.1) 

(4.2) 
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Equivalently, u G i?g(n) =: V, equipped with the norm | • |//i(o), satisfies 

a\iu,v) = {f,v) Vv G (4.3) 

where the bilinear form a\ : x i?o(n) —>■ K and the linear form / : i^o(^2) —> K are 


defined as 


ax{u,v) := / {ixu'v'+ auv) dx , {f,v)= / {fiv — f2v')dx. 

Jn Jn 


IQ JQ 

In view of the approximation of the operator A\ we introduce the metric space 

A:={X = {u,d)GL°°{n)xL^{n): < P < D*, - < d < a*} 

where b*, b*, S’*, ct* are positive constants defined as follows. Suppose that the pair {P,d) 
approximates {v, a) with error 




then it is easily seen that 


V* := — < V < V + — =-.v 


k-cr||io=(n) < y; 


— ^ ^ * 1 ^* —* 
-cr* — <a <a + y=:cr 


(4.4) 


Furthermore, using the Poincare inequality ||'c||^ 2 (q) < have 


{v* - ^cr*)k|^i(n) < ax{v,v) < (v* + 


ffi(n) 


for all V G Hq{Q), A G A. We conclude that setting a, := — ^cr* = and a* = 


V* + \a* = V* + \cr* + it holds 


V^\v\m(<rt) < III^^IIIa < v^l'cUqn) 'dv G VA G A 


(4.5) 


with |||f|||^ := ax{v,v). The space A will be a subset of A containing the coefficients A of 


the problem (2.1); it will be defined later on. 


Concerning the definition of the space F containing the right-hand side, we write 
/ = (/ij h) G iF'iyi.) X LA‘{VI) =: F (note that different couples in F may give rise to the 
same / G H~^{Vi)). 


We now discuss the /ip-discretization of (4.2). To this end, we specify that the binary 


master tree A is obtained from an initial partition, called the ‘root partition’, by applying 
successive dyadic subdivisions to all its elements. Later, cf. Property |4.H it will be needed 
to assume that this initial partition is sufficiently fine. Furthermore, with reference to the 
abstract notation of Section]^ given any (AT, d) G .ftx N we have p{d) = d. In consideration 
of this simple relation, throughout this section we will use the notation {K,p) instead of 
{K,d), i.e., the second parameter of the couple will identify a polynomial degree on the 
element K. We set 


Thus 


Vk,p = VpiK), Fk,p = Pp-i(if) X FpiK), 

Ak,p = {A = {P,d) G Pp+i(Ar) X Pp+i(iG) : < P < v*, — d^, < d < d*}. 

VS,={vG : v\Kn G ^PniKo) AD G 2)} 


will be the discretization space associated with the /ip-partition D. Furthermore, we have 
Fx> C F and Ad C A, with F and A defined above. The difference in polynomial degrees 
between the various components of the approximation spaces for data is motivated by 
the need of balancing the different terms entering in the local error estimators, see (4.18) 
below. 

At this point, we have all the ingredients that determine a Galerkin approximation as 
in ([2T^. 
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4.2 Computable a posteriori error estimator 


Given data (/d,Ad) G F® x Ad, let md(/d,Ad) G be the solution of the Galerkin 
problem ( |2.19 ) with such data. To it, we associate the residual r = r{uii, fn, Xjd) G 
defined by 


{r,v) = {f-r,,v) - ax^{uD{f-r,,\i,),v) Wv £ , (4.6) 


and satisfying (r, ud) = 0 for all vd G V,^. The dual norm of the residual is a natural a 
posteriori error estimator, since one has 

Ikliff-i(e) < llk(/D, Ad) - md(/d, Ad)|||ai, < lk||D~i(n) ; (4.7) 

in one dimension, such norm can be expressed in terms of independent contributions 
coming from the elements Kjj oi the partition 2), which are easily and exactly computable 
if, e.g., the residual is locally polynomial. To see this, let us introduce the subspace of 
HQ{n) of the piecewise linear functions on D, i.e., 

V^ = {v£ I e VF g 2)} C Kg 

and let us first notice that i7g(n) admits the orthogonal decomposition (with respect to 
the inner product associated with the norm |•|^rl(n)) 


© 0 H^{Kd) , 

neD 

where functions in Hq{Kd) are assumed to be extended by 0 outside the interval Ku; 
indeed, for any v gV, we have the orthogonal splitting 

V = VL+ ^ VKn , 

DeD 

where is the piecewise linear interpolant of on 23 and vko = ~ vl)\Kd ^ 

Hq{Kd). Recalling that {r,VL) = 0 for all vl G it is easily seen that the following 
expression holds; 

Del) 

where tkd denotes the restriction of r to Hq{Kd)- 

The computability of the terms on the right-hand side is assured by the following 
representation: for any 13 G 23, one has 

\\''^Kn\\H-){Kn) ~ \^Kd\'h){Kd)) 
where zkd € satisfies 

{z'Kn^^')L^(Kn) = (rKn^v) ^V£HI{Kd). (4.9) 

Writing ud = md(/d,Ad) and Ko = {a,b), and noting that, since / 2 ,d is a polynomial 
in Kd, 

{rKD^v) = {fl,D + f 2 ,D + - (Ti,UD)vdx = {rKD,v)L 2 {Ku)y (4-10) 

JKd 

it is easily seen that the solution zkd bas the following analytic expression 

ZKnix)= G{x,y)rKr,iy)dy , (4.11) 

JKd 
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(a-x){b-y) ^ ^ 

b-a ^ y 


where G'(a:, y) := < ~ is the Green’s function of our local problem (|4.9 1. 

[ b^ a; > y 

Thus, the squared norm Wi^KoWH-^iKn) local residual can be explicitly computed, 

since rx,:, is a polynomial. 

Summarizing, defining for any D G T> the local error estimator 


:= \zKn\m{KD) 


and defining the global error estimator as in (2.20), we have by (4.7) 

—= £d(ud(/d, Ad),/d, Ad) < |||u(/d, Ad) — ud(/d, Ad)|||ai 


< 


1 


/a* 


£d(ud(/d, Ad), /d, Ad), 


(4.12) 


(4.13) 


which in particular implies the reliability assumption (2.21). 


4.3 The module REFINE 

Hereafter, we present a realization of the module REFINE, that guarantees the discrete 


efficiency property (2.22), hence the contraction property of REDUCE. For every D G 


M C D the module raises the local polynomial degree to some higher value, whereas for 
D G D\M the local polynomial degree remains unchanged. No ^.-refinement is performed. 

To be precise, consider an element D = {Kx,pd) G M. Suppose that the local 
polynomial degree of the data is related to some p_D, in the sense that 

fl,D\Kn G IPpD-1 {Kd), f2,D\Kn G ^Pni^o), *^d|kd, O'dIko + 

Recall that md = md(/d,Ad) satisfies ud\ko G Fp^{K£)). Then it is easily seen that 
the residual r = r{uD, /d, Xd) is such that its restriction rx^ to Kx is a polynomial of 


degree px + Px + 1, while the function zx^, defined in (4.9) is a polynomial of degree 

Px-=PD+PD + i- (4.14) 

Therefore, the module REFINE builds T) = 2)(M) G = D with ©(M) > 2) as follows: 

{Kx,px) for U G M 


T) = {D} with D = 


D 


for U G CD \ M. 


In order to prove (2.22), consider a marked element D G M. Setting Fp^{Kx) '■= 
^PoiKx) F\ Hq{Kx) and recalling that zxd G Fp^{Kx) we have 


yD,D(wD,/ d. Ad) = \ZKn\H^(Ko) = ®Op 




{rKu,w) 




mm(XD) 


= sup 

RfGPO \W\m(KD) 


(4.15) 

On the other hand, the Galerkin solution = u,x{fx, Xx) is such that its resid¬ 
ual f = r{ux, Jt), Xx) satisfies {fxj^,w) = 0 for all w G ¥p^{Kx)- Thus, denoting 
by axj^^Xoi'G) the restriction of the form to H^{Kx) x H^{Kx), and setting 


1\It,,Ku = we get 

?7d.d(md,/d. Ad) = 


sup 


{rxu - rKo , w) 




sup 


aAi,,iA £,(md — ud,^ riTiii m 

-j—j- < Va |uD - Ux\1\t,,Kd ■ 
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Squaring and summing-up over all D S M, we obtain 


md,/d, Ad) < a*|||uD - ud 


llL 


(4.16) 


which immediately implies ( 2 . 22 ). 


Remark 4.1. The choice of the error estimator and the refinement strategy indicated 


above guarantees that the reliability assumption (2.211 and the efficiency assumption (2.22) 


are fulfilled, hence the eonclusions of Proposition 2.1 


hold true. Actually, one can be more 
precise, since using (4.13) and (4.16) and following the steps of the proof of Proposition 
\2.1[ we get that the sequence of Galerkin approximations built by a call of REDUCE 


satisfies the contraction property (2.25) with contraction factor k = y^l — 


4.4 Convergence and optimality properties of hp-AFEM 

In this section we discuss the convergence and optimality properties of our adaptive al¬ 
gorithm hp-AFEM in the present one-dimensional setting. To this end, we first specify 
the abstract functional framework introduced in Sect. We already set V := Ro(fl) and 
F := Lf{Tl) X Lf{Vl). Concerning the space A containing the coefficients of the operator, 
we assume stronger regularity than just L°°{Q) in order to guarantee that the piecewise 
polynomial approximations of the coefficients still define a coercive variational problem. 
To be precise, from now on we assume that A = {v, a) belongs to the space 


A := {A = (t) £ P[^{Q) X iJ^(n) : < v < v*, 0 < cr < cr*}. 


(4.17) 


Here, in view of (2.2), we choose to work with a smoothness space of Sobolev type with 


summability index 2 , so that squared best approximation errors are non-increasing under 
/i-refinements. We notice that it would be sufficient to require the coefficients to be 
piecewise on the initial partition. We decide to work under stronger assumptions just 
for the sake of simplicity. 

We now define the projectors Qk,p introduced in Sect. 


2.1 


To this end, 


let p 


L{L'^{K),Fp{K)) be the L^-orthogonal projection and n)^p £ L{H^{K),¥p{K)) be the 


iJ^-type orthogonal projection defined as follows: if u £ H^{K) with K = [a, b] then 


(a:) :=c+ f (H^^ (t) dt 
J a 


where the constant c is such that dx = dx. 

Then we define Qk,p £ xFx A,Pp(A) x (Pp_i(Ar) xPp(A)) x (Pp+i)/!") xPp+i(A)) 
by setting 


QK,p{v,f,\) := {Ii]^pV\K, ^K,p-lfl\K, ^K,pf2\K, ^K,p+lt^\K, ^K,p+iO'\k)- 
At last, we define the local error functionals eK,p. We set 

^K.pfv, f,\) := 1(1 — p)u|if -|- 5 ^osc|'p(/. A) 
where i5 > 0 is a positive penalization parameter to be chosen later and 

osc|-p(/. A) := II “(I ~ n5f_p_i)/i|_R-|li2(K) -I- 11(1 — n5(-_p)/2|A:|||2(^) 

+ 1(1 ~ ^]i,p+i)^\K\^m(K) + 1(1 ~ ^]<:,p+i)'^\k\’hi{k) 

where h = |i£|. Note that the choice of polynomial degrees is such that for smooth data 
the four addends above scale in the same way with respect to the parameters h and p. 


(4.18) 


(4.19) 
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Furthermore, the data oscillation that appears in (4.181 is of higher order with respect to 
the projection error for the function v. 


It is straightforward to check the validity of (2.2). We recall that given a partition 


2) S D, we denote by /d = (/i,d,/ 2 ,d) and Ad = (vvjO'd) the piecewise polynomial 
function obtained by projecting / and A, respectively, element by element as indicated 
above. Note that while /© G Fd C F, Ad need not belong to A. Given a partition D G D, 
we will set 

osc|,(/,A) := osc|,(/,A), 

dgd 

where osc|,(/, A) = osc^^ p^(/. A). 

The following result provides a uniform bound on the approximation error of the 
coefficients of the operator, assuring that Ad G A. 

Property 4.1. Let D be the root partition with polynomial degree equal to one on each 
element domain. Assume that 3C(T)) is sufficiently fine for the given data X € A, in the 
sense that for each K G 3C(CD) it holds 




Then for any B G D we have (4.4), i.e., 


W - i^D||L~(n) < y, 




\a - crD|U~(n) < y ■ 


Consequently, Ad G Ad C A. 

Proof. For any D = {Ko,Pd), let K G %{T)) the element of the root partition containing 
K]j. Then, we have 

\i^ - ^Kn,Pn + l')'^\Kn\HHKo) < \ 0 -- H^K) ^ y- 

On the other hand, set ip = {I — po+i)^li<'D! recalling that ip has zero mean-value 
in Ko, it vanishes at some point xq G since it is a continuous function. Writing 
ip{x) = ip{xo) + ip'{t)dt for any x G Kjj yields 

\lpix)\ < \x - Xo\'^^'^\\iP'\\l2(Kd) < I 

whence the result immediately follows after observing \Kd\ <1 . □ 


We now focus on the abstract assumptions (2.7)-(^2.9F 


Proposition 4.1. In the present setting, assumptions (2.8|-(2.9| hold true. Furthermore 


if 6 is chosen sufficiently small, then (2.7) is fulfilled. 


Proof. We start by verifying condition (2.9). For any v,w G Hq{XI) and for any B G 
and any B G B, it holds that 




Ko.PoJ^^IKdIH^Kd) 


= inf \w\Kn 

< inf 

¥>GPpQ (Kd) 


- PlmiKo) 

\v\Kd - TIhPKd) + l (^ “ w)\Ku\hPKd) 

= 1(1 - ^Kd,Po)'^\Ko\hi^{Kd) + K'C - w)\Ko\m{KD)- 
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Two applications of a triangle inequality show that 


Ed(u,/, A )2 - Ei,('u;,/, A )2 


< 


VdgD ^ 

— ^|(I — + (5 ^osc|)(/,A)j ^ 

I 

I (I ~ ^D^^KolH^iKn)) 


\DeX> 


< E l(i-nl,h 


^olHi(Kn) 


< U — W 


y> 


i.e., (2.9) holds true with constant C 2 = 1- 


Let us now verify assumption (2.8). Note that u(/i 3 , Ad) is well defined since Ad S A. 
Setting for simplicity u = u{f,X) and u = m(/d)Ad)) it is straightforward to check that 
u — u satisfies for any v G V 

ax{u — u,v) = {f — fD,v) — / {v — v-]x)u'v' dx — / {a — ax>)uv dx (4.20) 

Jo. Ja 

whence, using the Poincare inequality ||u||D 2 (f 2 ) < 2“5and selecting v = u — u, 
we obtain 

at.\u — u\h^{q) < ll/i —/i,D||i/-i(n) + 11/2 —/ 2 ,d||d 2 (q) 

+ - J2D|Uoo(n) + - \\cr - crD||L=°(n)^ |M|//i(n)- (4-21) 

We now bound the quantity on the right hand side of ( 4.21[ ) in terms of osc|,(/, A). To 
this end, starting with the first term, we have for any v G Hq^^I) 


ifl - fl,T>,v)L2{Q) — ((I - ^Kn.vn-l)fl\KD,v)L^(KD) 

D^D 

D^D 

- E 11(1 “n?f^_p^_i)/i|if^||L2(^^)||(i-n5f^_p^_;^)uiiG„)||D2(^^) 

(4.22) 


dgd 


By the classical /ip-error estimate for the orthogonal L^-projection upon ¥p^{Kj:,) (see, 
e.g., [3i Corollary 3.12] ) we have \\{l - ^Ku,PD-MKn)\\L^Ku) < C^\v\H^Kn) for 
some constant C > 0. Thus, we get 


||/i-/i.D||i/-i(0)<C' E ll/(I-nL.P.-i)/i|ifJlW.) • (4-23) 

Vdgd ™ / 


Concerning the second term on the right hand side of (4.21), we simply write it as 


11/2 - /2,D||L2(n) - I 11(1 - n?f^_p^)/2|_ffDllz,2(;f^) 

VnGD / 


(4.24) 
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Coming to the third and fourth terms, we first observe that 


|M|/i-l(n) <—(22 ||/i_d||j;,2(q) + ||/2,d||l2(q)^ 

— ( 2“2 + ||/ 2 ||i 2 (Q)^ =: C'(/), 


(4.25) 


< — 
a* 


since fi^D, i = 1,2 is locally an L^-projection of fi. On the other hand, using the same 
argument as in the proof of Property |4.1| we get 


liy - VTi\\L°-{n) - max 11(1 - 

- “I'd “ ^Kd,Pd + i)^\Ko\h^Kd) 

- ( \i^ - ^Kn.Vn + l)^\KD\Hi(Kn)j ' 

\DeD / 


(4.26) 


A similar result holds for ||cr — crD||L=o(f 2 ). Substituting (4.23)-(4.26) into (4.21) and 
recalling (4.19) we get 


1 

~ + C + 1^ ^ ^2 osci))/, A)^ . 


Thus, setting (7 := ^ (iC'(/) + C + 1) and recalling (|4.18|), we conclude that 


(4.27) 


i{f,X)-u{f^,X^)\H^n)<C6 ^ Y,(^d{wJ,X)] = CS"^ E^{w, f, 


(4.28) 




for any w € Hglfi). This proves that (2.8) is fulfilled with Ci = CS^. Finally, choosing 
any 6 such that Ci < 6 we fulHll (2.7). □ 


We conclude that choosing 6 sufficiently small we may apply Theorem |2.1[ This 
leads to the conclusion that for solving (4.2), where / = (/i,/ 2 ) G L?{^) x T^(n), and 
A = (z 2 , (t) € A dehned in (4.17), and with a root partition 2) that is sufficiently Hne such 
that it satisfies Property 4.1 hp-AFEM is an instance optimal reducer, in the sense of 
Theorem 12.11 of the error functional 


Ad(w(/,A),/,A) = inf \u{f,X)\K^-V?m(Ko)+^ ^osc|,(/,A), 


over all D S D, where osc|j(/, A) is defined in (4.19). 

Finally, we consider assumption (2.14). At first, we note that in one dimension all 
partitions are trivially conforming, i.e., D'^ = D. Next, we observe that the following 
result holds. 

Lemma 4.1. For any 2) G D and any v G 77g(f2) there holds 

|2 _ • r I , ,„|2 


inf \v - W']:)\i,i(n) = 7 inf In 


- fIh^Ko)- 


(4.29) 
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Proof. For D gT), let qd G Fp^{Kd) be such that |w - golH^Ko) = inf<^GPpo(^ d) ^ “ 
p\h^{Kd)- Define g G L 2 {fl) by g\Kr, = q'd D gT>, and wd G H^{^) by w-r>{x) = 

g{s)ds. From v', we infer that wd(0) = wd( 1) = 0, and so w-d G Vig. 

Moreover, |n - = EdgD “ (Id^h^KuY ° 


Observing that 

, 1 ^ “ P?m{Ko) = 1(1 “ ^]iD.PD)'^\KD?m{KD) - ^D{v,f,X) 


for any / G F", A G A, we obtain the following result. 
Proposition 4.2. For all D G V) and all v G iFQ(O), one has 


inf \v-WD\HUn)< inf ED(n,/,A)2, 
KJBeFg (/,a)gfxA 


i.e., for C := / assumption (2.14) is fulfilled with = 1. 


As a consequence, (2.16) and (2.18) are fulfilled with Ca = (74 = 1. Since hp-AFEM 
calls the routine REDUCE with the hxed value g = , and by Proposition 

the number of iterations in REDUCE is bounded by O(logp“^), we are guaranteed 


2.1 


that the number of iterations performed by REDUCE at any call from hp-AFEM is 


uniformly bounded. On the other hand, recalling (4.14), for each iteration in REDUCE 


the polynomial degree in each marked element is increased by a constant value depending 
only on the local polynomial degree in the input partition. Thus, even in the worst-case 
scenario that at each iteration all elements are marked for enrichment, we conclude that 
the output partition of REDUCE has a cardinality which is bounded by a fixed multiple 
of the one of the input partition, which is optimal as it is produced by hp-NEARBEST. 

Another obvious, but relevant application of Lemma 4.1 is that hp-AFEM is an 


instance optimal reducer over (D G D of the error functional written in the more common 
form 


inf 

wx,&V£ 


i(/) ~ ^ 


OSCr 


,(/,A). 


5 The Poisson problem in two dimensions 

On a polygonal domain O C K^, we consider the Poisson problem 

f —Am = / in O, 

\ u = 0 on 90, 

in standard variational form. We consider right-hand sides / G Lf{Vl), and so take V = 
iFg(O), F = Lf{FL), and A = 0. We equip Ro(0) with | • and F[~^{Ft) with the 

corresponding dual norm. 

Let OCq be an initial conforming triangulation of O, and let in each triangle in Xq one 
of its vertices be selected as its newest vertex, in such a way that if an internal edge of 
the triangulation is opposite to the newest vertex of the triangle on one side of the edge, 
then it is also opposite to the newest vertex of the triangle on the other side. As shown 
in [3 Lemma 2.1], such an assignment of the newest vertices can always be made. 

Now let IK be the collection of all triangulations that can be constructed from Xq 
by newest vertex bisection, i.e., a repetition of bisections of triangles by connecting their 
newest vertex by the midpoint of the opposite edge. With each bisection, two new triangles 
are generated, being ‘children’ of the triangle that was just bisected, with their newest 
vertices being defined as the midpoint of the edge that has been cut. The set of all triangles 
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that can be produced in this way is naturally organized as a binary master tree having 
as roots the triangles from Xq. The triangles from ^ are uniformly shape regular. The 
collection K of triangulations of 17 is equal to the sets of leaves of all possible subtrees of 

For K G we set Vk = H^{K) and Fk = L‘^{K), and for d S N, we set 


VK,d ■= Vp{d)iK), Fka := '^p{d)-i{K), 


(5.1) 


with, as in Sect. 2.1 p = p{d) being the largest value in N such that dimPp_i(i4r) = 
example, for d = 1,..., 10, we have p = 1,1, 2, 2, 2,3, 3, 3, 3,4. 

Remark 5.1. Alternatively, one can select sequences of strictly nested spaces {Vx^dfdi 
{FK.d)d with the condition that for the values of d of the form for some p =: 

p{d) G N, definitions in (5.1) hold. 

For D = {Kd, do) G 3C x N, we write Vd = VkoAd^ = FkoAd and po = pido)- 
Note that with the current definition of Vd, this space is uniquely determined by specifying 
Kd and pd- For some constant d > 0 that will be determined later, we set the local error 
functional 

11/ ~ /£>l!i2(K^)j 


eD{w,f):=eD{w) + 5 inf 

Pd foev^n-AKn) 


where 


We define 


CDiw) := inf 

{uidGPp^CATd): w} 


\w -WD\m(Ku)- 


Qd{wJ) := {wD,fD) 


(5.2) 


(5.3) 


as the pair of functions for which the infima are attained. 

Having specified the master tree A, the local approximation spaces Vd and Fd, the 
error functional eD{w,f), and the projection QD(w,f) = {wD,fD), we have determined, 
according to Sect. |2.1[ the collection of dp-partitions D, the approximation spaces Vd and 
Fd for 2) G D, the global error functional 


Ed(w,/)= eD{w)+S ^osc|,(/), 

DeD 


(5.4) 


where 




as well as the projection /d := YIdgTi fo- 


We proceed with verifying assumptions (2.7), (2.8) and (2.9). 
Proposition 5.1. There holds 

sup |Ed(w,/)5 — Ed(u,/) 5| < lire — u||v V2) G D, \/v,wGV, 
/6P’ 


i.e., (2.9) is valid with C 2 = 1- 

Proof. For v,w G V, it holds that eD{w)^ < eD{v)^ + |u — which yields the 

proof using the same arguments as in the proof of Proposition 4.1 □ 
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Proposition 5.2. There holds 


Hf)-u{M\m^n)<VS inf EdK/)5 VB g D, \/f g L^{n), 

weH^id) 

i.e., (2.8) is valid with Ci ~ '/S, and, when S is chosen to he sufficiently small, so is (2.7). 

Proof. Since / i-G u{f) G T/q (il)) is an isomorphism, it is enough to estimate 

11/ ~ /D||H-i(n)- To this end, we note that for K being a triangle and p G N, it holds that 


m 


sup 


inf 


\w - ^ diam(it:) 


\w\h^{k) ^ p+1 

only dependent on a lower bound for the smallest angle in K. Consequently, we have that 

inf„g_Fj,(/ - /d,w - v)i2(f2) 


\\f - IdWh-PQ) = sup 

weH^in) 




< 


sup 


E 


DeD 


\f^D\ ^ 


11 / - fD\\L2{KD)\MH^KE 




< ^osc^(/). 


(5.5) 


□ 


5.1 Conforming /t-partitions, and conforming hp finite element 
spaces 

For the design of a routine REDUCE, in particular, for a posteriori error estimation, it 
is preferable to work with /i-partitions that are conforming. Let 


:= {3C G K: 3C is conforming}. 


As shown in [6l Lemma 2.5], for 3C G K, its smallest refinement G K'^ satisfies < 
#X. 

With the subclass 

:= {T> G D: 3C(D) G K^j, 

we define C : D —^ D'^ by setting 6(1)) = I), where 1) is defined as the partition in D° with 
minimal #1) for which 1) > D. That is, = X{T>Y, and pn = Pd for D € 2}, H S D 

with Kjj (Z Kjj. 


Unfortunately, sup^g 


#e(D) _ 

#D 


= 00 , I.e., 


(2.18) is not valid. Indeed, as an example, 
be such that Ki G X{1)), 


consider Xq to consist of two triangles Ki and K2. Let D G 
with corresponding polynomial degree p{d), and that in X{T)), K 2 has been replaced by 
2^ triangles of generation N, each with polynomial degree 1. Then ffT) ~ d + 2^. Since 
1C(6(D)) = X{'Dy contains in any case ~ 2^^"^ triangles inside Ki, so with polynomial 
degrees p{d), we conclude that #6(D) ^ 2^ + 2^/'^d. By taking say d ~ 2^, we conclude 
the above claim. 


The fact that (2.18) does not hold implies that, unlike for an /i-method, we will not 


have a proper control on the dimension of the finite element spaces that are created inside 

REDUCE. 


From (2.6), recall the definition Ug = Vb n Hq{Q) for T> G D'^, and from (5.2)-(5.3), 


recall the definition of eniw) and wd for D G B and w G H^{Ko)- The main task in this 
section will be the proof of the following result. 
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Theorem 5.1. Setting, for 2) S D, ||pd||oo := Pd > for T) it holds that 


inf < (1+ log|bi,||oo)^ ^ 


Since, for D G D, obviously J^D^eiD) ^d{w) < ^oiw), Theorem 5.1 implies 

C3,D ~ (1 + log IIpdIIoo)® ■ 


(2.14) with 


For an underlying /i-partition that is conforming, Theorem |5 .1 1 says that the error in 
il^-norm of the best conforming /ip-approximation of a w G iJQ(n), is at most slightly 
larger than the error in the broken il^-norm of the best nonconforming hp-approximation. 

The proof of this remarkable result will be based on Veeser’s proof in [45] of the 
corresponding result in the ‘/I’-setting. In [1^, the result is shown by taking v to be 
the Scott-Zhang ([ID]) quasi-interpolant of w. This Scott-Zhang quasi-interpolation is 
constructed in terms of the nodal basis, and the proof relies on an inverse inequality 
applied to these basis functions, which inequality involves a multiplicative factor that is 
known to degrade seriously, i.e. not logarithmically, with increasing polynomial degree. 

In our proof the role of the nodal basis on a triangle will be played by the union of 
the three linear nodal basis functions associated to the vertices, the polynomials on each 
edge that vanish at the endpoints, that will be boundedly extended to polynomials on 
the interior of the triangle, and, finally, the polynomials on the triangle that vanish at its 
boundary. We will construct a IId G £(7Jq (fl), Idg) such that, with 


Pd d '■= max po' VZ? G 2), 

’ {D’^DiKoinKoi^^i} 


(5.6) 


it holds that 

\iw-^Dw)\Kn\HHKu) ^ (1 + logPD.I))^ 


eD'{w) WDeTi, (5.7) 

{D'GD:ifD/nifo#0} 


which obviously implies the statement of the theorem. Since the right-hand side of (5.7) 
vanishes for w G Idg , because it even vanishes for w G Vd , the mapping IIii is a projector. 


Proof. {Theorem 5.1) Let D G D'^. In order to show ( |5.7[ ), it is sufficient to show 

\{TI'dw)\kd-wd\hi(Kd) - (l+log7'£’.D)^ E eo'iw) VL> G 2), w G i7o(n). 


{D'e'D:Kj^,nKo^ll} 

(5.8) 

Let 3M(2)) and £(2)) denote the collection of vertices (or nodes), and (closed) edges of 
X{T>). To construct IId, for e G £(2)) we set 

Pe,D '■= min{p£i : £> G 2), e C BKd}, pe^D := max{p£) : D G 2), e C dKo}. (5.9) 

With the mesh skeleton dX{T>) := UDGDdKjj, we set 

Vav ■■= {v G C(52C(D)) : v\e G Pp,,^(e) Ve G £(2))}, 

VaD ■■= {V G C(dX(V)) : vU G Pp,.^(e) Ve G £(2))}. 

We construct IlaD G Zl(neg£(ii) H^{e), Van), and an auxiliary Uav G '^dleGECD) H^e),VaD), 
such that 


(naDw)| 9 n = 0 for all v = (ue)eg£(D) G H^^{e) with z;e = 0 for e C 912. (5.10) 

eG£(D) 
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For any triangle K with edges ei, 62 , 63 , there exists an extension Ej^ G {OK), H^{K)) 

that, for any p G N, maps C{dK) n ni=i IPp(ei) into ^p{K) (see e.g. [21 Sect. 7]). Defining 
flu by 

{UDw)\Kr, ■= EKo{(^dT)W\dx(D))\dKD) +'>^D - EKriiwoldKo)^ (S-H) 
in view of the definition of and ( |5.10[ ), we have IId G tdg). 

To construct IlaD, tid'D, for each G 7^(2)) we select some 

G £(B) with V & e,j and C dVl when v G d^. (5-12) 

For V G 7^(25), by (f)^, we denote the nodal hat function, i.e., (f)^, is continuous piecewise 
linear w.r.t. 2C(2)) and 4>v{v) = ^v,v G 2^(2)). For e G £(23), let 

Oe : i25(e) —>■ H^{e) be the (e)-orthogonal projector onto Pp^ ^, (e), 

Qo.e : 22 ® (e) —>■ 22 5 (e) be the 225 (e)-orthogonal projector onto Pp^ ^, (e) n 22 q (e), 

Qo.e : 222 (e) —>• iJ 2 (e) be the 222 (e)-orthogonal projector onto Pp^ ^, (e) n 22 g(e). 

Denoting the endpoints of an e G £(23) by vi^e,V 2 ^e, we now define UdD and by 
setting, for v = {ve)eeE{D) e UeeEii}) 


Z Z 

(ffp!D^) |e ■“ ^ P ^ 

i—1 i—1 

2 2 

(n^2)^) |e ■ ^ ^ (^Qeu^ P ^0,e^^e ^ ^ 


for any e G £(23). It is clear that Tid'o maps into VgD, and, thanks to (5.12), that it 
satisfies (5.10). Similarly, nan maps into VgD 

These definitions show that, for D G 23, {I\.dw)\kd depends only on w|u{j^^,.£ 3 /gD, if£,/ni<'D/0}- 
Therefore, in order to prove (5.8), a homogeneity argument shows that we may assume 
that Kjg is a uniformly shape regular triangle with 


|22r,| = l. 


Since the extension 23i^Q G L{H2 (dKo), H^{Kd)) can be chosen to be uniformly bounded 
over all such Kjg, in view of (5.11) in order to arrive at (5.8), and so at the statement of 
the theorem, what remains to be proven is that 

\\{^gDw\d'X{'D))\dKD 

< (1 + \0gPD,Df 


'H2{aKD) 

e_D'(w) VD G 23, lu G 22^(22). 


(5.13) 


{D'(ED:Ko,nKr,¥^<l)} 


In [21 Thms. 6.2 and 6.5], it was shown that on an interval 2 of length ~ 1, it holds 


that 


lkllLoo(/) ^ (1 PlogP)®ll^lljji(P) V0GPp(2), (5.14) 

ll^ll 1 <(l + logp)||z||^i VzGPp(2)n22i(2). (5.15) 

These estimates will be used hereafter. 


Lemma 5.1. For v G 22(23) n dKjg, e, e' G £(23) with = v, we have 

\{QeW\e - Qe'W\e>)iiy)\'^ < {I + log PD')eD'M VwG 22^12). 

{D'GB: 
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^n—1 



Figure 2: Notations relative to the proof of Lemma 5.1 


Proof. Consider the notations as in Figurej^ Using that for 1 < i < n, {Qf,.WDi\ei)ii') = 
(Qei_iWDjei_i)(j^), we have 


((5e„w|e„ - QeoW\eo){'^) = 

n—1 

{pe^iw -WKoJ\e„ + ^ Q eiiw K - W + W - WKo,)\ei + Oeo Uo) (^)> 


and so, using (5.14) and the trace inequality, 


|(Qe„w|e„ - QeoW\e„){v)\ [(1 + logPe._i ) " + (1 + logPe J " ] 6^. (w) ^ . 

2=1 


We continue with the proof of Theorem |5.1| As a first application of this lemma, 
we show that it suffices to prove (5.131 with Hg-r, reading as tlg'T). To this end, for 
e G £(®) n dKjy, let U' G 2) such that e C dKiji and Pe,i) = Pd'- Then 

2 

(naDW^|a3C(D))|e — (naDUi|a3c(D))le = {Qo,e — Qo,e)(fw\e — 

i=l 

2 

— {Co,e Qo,e) I e '^D'\e ^ ^ ' 
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From (5.15), the trace theorem and the property || i, < 1, we infer that 


\{^aT)W\d‘X{D)\e—{^dDw\Q-x{^)\e 


Hl(e) 


A 

< (1 + \0gpe,D)(eD'{w)^ + ^ “ 'ICD'|e) (t'i.e) |) ■ 


2 = 1 


(5.16) 


Writing 


Qe^ w;|e„. -WD'\e = Qe^. w\e^. - Q eMe + Q e{w\e - W D>\e), 


and applying (5.14) as well as the trace theorem, shows that 

<\[Qe^. W\e,, - QeW|e)(j^j,e)| + (1 + l0g(Pe,D))5eD'(lc)5- 


(5.17) 


By combining (5.16) and (5.17), and applying Lemma |5.1| to the first term on the right- 
hand side of ( |5.17| ), we conclude that 

II {{^d'D-^d'D)w\dOC(V))\dKD 

< (1-I-logp_D,D)^ ^ eD'{w) VD G T), w G 

{D'Gl':ifj3/niCD/0} 


(5.18) 


As a consequence, what remains to show is (5.13) with Hqt) reading as llaD, that is, 
to show that 


||(naDu;|a3c(D))|aifi, - 

< {1 + \ogPD,T>)^ ^ CD'iw) VD G D, w G 

{D'GD:iCD/nifD 5 ^ 0 } 


(5.19) 


Let us first consider the situation that C dKjj for all ^ G 7^(2)) n dKo- Then 
((/ — fidi>)wD\docCD))\dKD — 0 (this is generally not true for flan), and so 

\\{fldT>w\Qx(V))\aKn ~ = Wi^dDiw - 'WD)\aKn)\dKD\\ (5-20) 


To bound the right-hand side, let us write v = {w — WD)\dKo- I^or edges ei, 62 of dKu, 

(5.21) 


and := Cl n 62 , an application of (5.14) shows that 




For an edge e C applications of (5.15) and (5.14) show that 

2 




II Oo,e ^"^le ^ i ^ |e^ I 

2 = 1 

2 

< (1 -f logPe,D)||(5o,e('y|e “ '^{Qe^^ .v\e^^ J{iye,i)(/}u,Je'j \ 

2=1 

< (1 -f logPe,D) (||?^|e||j^i + max(l -h logPe,.^ .,d) 5 ||u|e„^ 


(5.22) 


H2(e) 




Combination of (5.20), (5.21), and (5.22), together with an application of the trace theo¬ 
rem, show that, in the situation of Cj, C DKd for all v G 14(1)) n dK^, 

\\{ndTiw\dKo)\dKo -WDlaKoWlritfi,. , ^ (1 + logp_D,D)^e_D(w) G H^{Kd), 

H2 (dKn) 
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which implies (5.19). 

Consider now the situation that for one (or similarly more) v G (11(1)) H dK^, e^, <f. 
dKi). We estimate the difference, in i /2 (^iiT^))-norm, with the situation that is equal 
to some edge e C dKo. Applications of ( 5.15| ) and Lemma 5.1 show that 

II XI {I-Qo,e){{QeMe.-QeW\e){v)(l}^\e)\\^Xi^g^^^ 

{eG£(D)n9ifD:e3j/} 

^(1+ X {^ + '^Ogpe,-D))\{Qe^w\e„-QeW\e){v)\ 

{eG£('D)n3/C£j:e9i>'} 

< (1 + logPD,D)5 X eD/(w)5, 

{D'GD:iC^,niCD#0} 


which completes the proof of (5.19), and thus of the theorem. 


5.2 The routine REDUCE 


For 2) G D'^, with wd we will denote the best approximation to w from = Ed C Hq{VI) 
w.r.t. I • |i/i(n)- For w = u{f), being the solution of the Poisson problem with right-hand 
side /, ud(/) turns out to be the Galerkin approximation to u{f) from E^. 

In this section, we will apply results from |l{5j on residual based a posteriori error 
estimators in the hp context. These results were derived under the condition that the 
polynomial degrees po and for D,D' G D G W with Kd H ^ 0 differ not more 
than an arbitrary, but constant factor. Fixing such a factor, let D'^ denote the subset 
of those 2) G D“ that satisfy this condition. Obviously, for each 2) G D“, there exists a 
2) G D° with 3C(f)) = 2C(2)) and (D > 2). Unfortunately, even for the smallest possible of 
such 2), let us write it as 22(2)), the ratio cannot be bounded uniformly in 

2) G DC 

In view of the replacement of by D'^, the mapping C : D —)■ D'^ constructed in the 
previous subsection has to be replaced by C := D —>■ D° : 2) 2)(C(2))). From now on, we 

will denote D'^ as D'=, and C as C. Since obviously 2) can be constructed such that \\p^ ||oo = 

IIpdIIc 


(1 -b log(llpDlloo))^ 


with these new definitions (2.14) is still valid with Cs^d 
and, as before, unfortunately supDgn #^(F’)/#® = oo- 

We note that in the present application, for 2) G D°, /d G Fd, and a desired reduction 
factor g G (0,1], REDUCE(p, 2), /d) has to produce a 2) < 2) G D'^ such that |u(/d) — 
^©(/©(Iffqn) ^ 0 |w(/d) ~ WD(/D)|//i(n)- As explained in Section 2.3 the Eth iteration 
of hp-AEEM performs a call of REDUCE( ^ , 6(2)^), /dJ. The scalars /i 

and w are parameters as set in hp-AFEM. They depend on the constant b from hp- 
NEARBEST, cf. Sect|2.2[ the constant C 2 , here bein g eq ual to 1, cf. Proposition |5.1[ 


and the constant Ci, here being ~ \/S, see Proposition 


5.2 


The scalar ^ is a parameter 
in the definition of the error functional E, see ( |5.4[ ), that is chosen such that C 1 C 2 < b, 
cf. (2.7). The only possible dependence of the required reduction factor ^ on 

Di is via the value of As we have seen, Ca^Di ~ (1 + log ||pDi||oo)E meaning that 

when the maximum polynomial degree in 2), tends to infinity, this reduction factor tends 
to zero, but only very slowly. 

The construction of REDUCE will follow the general template given in Sect. |2.4| 
We will verify the assumptions (2.26), (2.27), and (2.28). For (w,/) G x 


2) G D'^, and Z) G 2), we set the residual based (squared) a posteriori error indicator 


’7i),d('^)/) ^ 2 + ^'^\\‘l^{Kd) + 

Pd 


X 


{eeE{'D)-.eCdKDnn} 


2pe 


-||[Vw • nelll 


2 

L2(e)) 
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where Pe,D is from (5.9), and define 


£ d ( w ,/ d ) := ( XI ^D,viwJv)j 

VDeD / 


1/2 


The following theorem stems from [551 Theorem 3.6]. Inspection of the proof given 
therein shows that the local lower bound provided by the (squared) a posteriori error 
indicator applies to any w G and so not only to the Galerkin solution. 

Theorem 5.2 (‘reliability and efhciency’). There exists a constant i? > 0 such that for 
D gW and /d € Fd, 

|u(/d) - MD(/D)|ffi(n) < ^£i)(md(/d), /d). (5.23) 

IIpd such that for all 


For any e > 0, and all D G there exists an rn^ 
fx) G Fx), and w G Vfy, 


rx),e£viwjx) < |m(/d) - w^lffi(n)- 


(5.24) 


Corollary 5.1 (‘stability’). With rx>^e 'In Theorem 5.2 for all D G D°, fx G Fx, and 
v,w G Idp, it holds that 




£a3 (^',/d) - £d(w,/d) 


< \v — w\ 


ffi(n)- 


(5.25) 


Proof. A repeated application of the triangle inequality, first in £2 sequence spaces and 
then in function spaces, shows that 

|£d(^^,/d) - £d(w,/d)| 

\DeX> {eG^(D):eC^ifDnn} 

_ 1 


where the last inequality follows from an application of (5.24) with “/d” reading as 0, 
and “ic” reading as v — w. □ 

What is left is to establish the ‘estimator reduction by refinement i.e. ( |2.28| ). Given 
M C D G D'^, we define !D(M) G D'’ as follows: The h-partition DC()f)(M)) is the smallest 
in K'’ in which each Kx for I? G M has been replaced by its four grandchildren in .ft; and 
for D G I)(M), it holds that po = Pd' where D' G T) is such that Kx' be either equal to 
Kx, or its ancestor in 3C(2)). 

Proposition 5.3 (‘estimator reduction by refinement’). For M C B G D'^, and I)(M) G 
defined above, it holds that ^B(M) ;< ffFl. For any fx G Fx, the estimator reduction 
property (2.28) is valid for 7 = 5 . 

Proof. This follows easily from the fact that each Kx {D G M) is subdivided into four 
subtriangles that have equal area, that each e G £(M) is cut into two equal parts, and 
that the jump of the normal derivative of rc G Vf, over a newly created edge, i.e., an edge 
interior to a Kx for D G B, is zero. □ 


Given B G O'” and / G Fx, let B = Bq < Bi < 


C 


be the sequence of 


hp-partitions produced by REDUCE(/), B,/d). We have established (5.23), (5.24), and 
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(5.25), for any fixed e > 0, as well as Proposition 5.3 Observing that |bi)(M)l|oo = Ibull 


an application of Proposition |2.2| now shows that in each iteration the quantity 


b(/D) - MDi(/D)|^i(n) + (1 - V7)»'D,e £d, ) /d)) 


where 7 = (1 — d)+ 0 / 2 , is reduced by at least a factor 1 — and that this 

quantity is equivalent to b(/D) — wi)b/i>)lffi(a)- view of ~ lbi>ll«D^~^^) 
conclude that in order to reduce the initial error b(/D) — UD(/i))bi(o) by a factor g by 
an application of REDUCE, the number of iterations that are required is 


M^log(l/p)|bD||L+^^ 


Remark 5.2. This result is not satisfactory because the number of iterations grows more 
than quadratically with the maximal polynomial degree. Yet, it improves upon the result 
stated in m, where the number of iterations scales with the fifth power of the maximal 
polynomial degree. 
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